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We study the growth of subhorizon perturbations in brane-induced gravity using perturbation 
theory. We solve for the linear evolution of perturbations taking advantage of the symmetry under 
gauge transformations along the extra-dimension to decouple the bulk equations in the quasistatic 
approximation, which we argue may be a better approximation at large scales than thought before. 
We then study the nonlinearities in the bulk and brane equations, concentrating on the workings of 
the Vainshtein mechanism by which the theory becomes general relativity (GR) at small scales. We 
show that at the level of the power spectrum, to a good approximation, the effect of nonlinearities 
in the modified gravity sector may be absorbed into a renormalization of the gravitational constant. 
Since the relation between the lensing potential and density perturbations is entirely unaffected by 
the extra physics in these theories, the modified gravity can be described in this approximation 
by a single function, an effective gravitational constant for nonrelativistic motion that depends on 
space and time. We develop a resummation scheme to calculate it, and provide predictions for the 
nonlinear power spectrum. At the level of the large-scale bispectrum, the leading order corrections 
are obtained by standard perturbation theory techniques, and show that the suppression of the 
brane-bending mode leads to characteristic signatures in the non-Gaussianity generated by gravity, 
generic to models that become GR at small scales through second-derivative interactions. We 
compare the predictions in this work to numerical simulations in a companion paper. 



I. INTRODUCTION 

The understanding of the current acceleration of 
cosmic expansion is one of the most pressing chal- 
lenges in cosmology. While the cause of acceleration 
is presently not well understood, a substantial num- 
ber of observations are being planned in the near 
future that will try to solve this puzzle. 

Cosmic acceleration can be explained in at least 
two different ways. Since in general relativity (GR) 
gravitational effects can be generated by pressure 
in addition to energy density (being, respectively, 
space-space and time-time components of the stress- 
energy tensor), the so-called dark energy solution to 
cosmic acceleration postulates a new contribution 
to the stress-energy tensor that dominates at late 
times and effectively has an equation of state where 
pressure is comparable to minus energy density (in 
natural units), leading to a repulsive force that fu- 
els acceleration. The simplest, best motivated, such 
model in agreement with current observations is the 
cosmological constant. 

Another possible explanation is that cosmic accel- 
eration signals a breakdown of GR at present cos- 



mological scales, rather than a new contribution to 
stress-energy. In this case, gravity is modified in 
the infrared (large-distance) regime to provide ac- 
celeration even though the universe is dominated by 
non-relativistic matter at late times. 

Mapping the expansion history of the universe 
from deceleration to present acceleration typically 
cannot distinguish between these two possibilities, 
as e.g. the Friedmann equation in a theory of mod- 
ified gravity can always be well approximated by a 
GR Friedmann equation with a dark energy com- 
ponent with suitable equation of state. The key is 
to add information on the growth of structure [1], 
which results from a competition between gravita- 
tional attraction and the expansion history: different 
theories of gravity will result in different growth of 
structure for a fixed (observed) expansion history. 
By now, a large number of tests along these lines 
have been proposed in the literature [1-14, 14-26], 
and will be pursued in upcoming surveys. 

Modifying gravity at large scales is a difficult task, 
and so far there has been limited success in finding 
consistent theories. GR is a very constrained theory, 
since it follows from requirements of locality in 3-1-1 
dimensions, massless spin two graviton and general 



covariance [27-31]. Any modification to it entails 
new degrees of freedom. When allowing for com- 
plex models of dark energy, however, the meaning 
of what constitutes a modification of gravity can 
become nontrivial to define. For example, many 
theories usually considered in the cosmology liter- 
ature as modification of gravity can be written as 
GR plus a (possibly non-minimally coupled) scalar 
field, i.e. some form of dark energy (a notable ex- 
ample in this class are f{R) theories [32, 33]). In 
such scalar-tensor theories, where only a scalar de- 
gree of freedom is added to GR, it is a matter of 
choice how to describe the physics (since both are 
equivalent). A more restrictive definition of grav- 
ity modification is to define departures from GR as 
such theories in which tensor perturbations propa- 
gate differently than in GR even in the absence of 
any contributions to the stress energy of the back- 
ground, i.e. gravitons behave differently than in GR 
even in flat space [34]. In this case one considers as 
gravity modification the addition of scalars only if 
they change the properties of the spin-2 sector, e.g. 
altering the speed of propagation of gravitons. 

Finally, the most restrictive definition of modi- 
fied gravity (and the one we adopt in this paper) in- 
volves modification in the number of degrees of free- 
dom in the spin-2 sector as a requisite, e.g. extra 
graviton polarizations, or many gravitons. Exam- 
ples include theories that realize the idea of massive 
gravity, e.g. large-scale extra-dimensions [35-38], 
degravitation [39] and bigravity [40, 41]. These the- 
ories generically have a scalar-tensor regime, where 
the predictions can be approximated by the addition 
of one or more effective scalars (e.g. the zero-helicity 
component of the graviton/s), but the theory in full 
is not equivalent to an arbitrary dark energy model. 
The best-known example is the DGP model [35], in 
which the large distance behavior involves from the 
4D point of view an infinite number of massive spin- 
2 states that describe the fact that the theory is 
fundamentally five-dimensional. This class of theo- 
ries also have more practical advantages. As already 
mentioned, general covariance imposes very power- 
ful constraints on the nonlinear couplings of spin-2 
fields. Requiring consistency leads to strongly con- 
strained theories that are easiest to test (or rule out), 
since they contain only a small number of free pa- 
rameters (instead of free functions). See e.g. [42] for 
more discussion along these lines. 



From the point of view of cosmology, the extra de- 
grees of freedom in the metric that behave as scalars 
(in the sense of the 3-l-l-SVT decomposition, see 
below) are the ones that lead to readily observable 
consequences. To be a viable theory of gravity, how- 
ever, any extra scalar degree of freedom must be 
suppressed in the solar system, where classical tests 
(light deflection in combination with non-relativistic 
motion) have shown GR to be a very good approxi- 
mation. 

To date perhaps the most promising mechanism 
for such a suppression is a generalization of the Vain- 
shtein effect found in massive gravity [43] , where ex- 
tra degrees of freedom become nonlinear through 
their self- interactions and decouple from matter, 
which then leads to the recovery of GR [123]. The 
Vainshtein mechanism happens under rather generic 
conditions and plays a role in various theories of 
modified gravity [34, 44-46]. Being a nonlinear ef- 
fect, the scale at which it happens is related to the 
amplitude of the sources; for stars such as the sun 
the Vainshtein scale is tiny by cosmological stan- 
dards, but for cosmological sources it is expected to 
be in the interesting few Mpc range [47] . Therefore, 
at scales smaller than the Vainshtein scale these the- 
ories will be consistent with GR, and thus difficult to 
distinguish from dark energy by cosmological obser- 
vations, i.e. the expansion history plus growth tests 
mentioned above would fail [124]. 

Such considerations show that probes of cosmic 
acceleration that are sensitive to small scale pertur- 
bations may yield no (or weak) tests on the cause 
of the expansion history, and could be mistakenly 
interpreted in favor of dark energy if the possibil- 
ity of nonlinear effects in the modified gravity sec- 
tor is ignored. It is therefore important to look for 
signatures of such nonlinearities in modified gravity 
theories that can be used as a diagnostic in future 
observations. That's in part the motivation for this 
paper and its companion paper II [48]. 

In this paper we study the growth of structure in 
brane- induced gravity in five dimensions [35] , known 
as the DGP model, which presents all the ingre- 
dients discussed above. This theory belongs to a 
broad class of modified gravity theories due to large- 
scale extra-dimensions [36-38], and presents simi- 
larities with other realizations of the idea of mas- 
sive gravity such as degravitation [39] and bigrav- 
ity [40, 41]. For theories that implement massive 
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gravity through something akin to the Higgs mech- 
anism see e.g. [49, 50]. 

The DGP model has two branches of cosmological 
solutions, the so-called normal and self-accelerating 
branches. The latter, in the presence of non- 
relativistic matter, makes a transition from early 
time deceleration to acceleration [51, 52] with de Sit- 
ter expansion in the long-time limit. In this paper, 
we concentrate on the self-acelerating branch as an 
example of a theory of modified gravity that leads 
to acceleration, even though several lines of work 
have shown that this branch has significant prob- 
lems related to the appearance of a ghost mode at 
the linear level [34, 53-58], and thus it may not be 
regarded as a viable description of our universe at 
all scales, particularly when quantum effects are in- 
cluded [125]. It's worth pointing out, however, that 
a ghost is not a necessary ingredient to obtain self- 
accelerated solutions [34]. 

On the other hand, as we discuss below, such is- 
sues do not seem to impact the cosmological solu- 
tions we study here, in the sense that cosmological 
perturbations that respond to a standard fluctua- 
tion spectrum source do not show any pathologi- 
cal behavior. In addition, and most importantly, 
it is expected that the techniques developed here 
will be widely applicable to the study of cosmologi- 
cal perturbations in other (consistent) large-distance 
modifications of gravity. We thus take the DGP 
model in the spirit of the simplest example of the- 
ories (parametrized by a single free parameter, as 
ACDM dark energy) where observational signatures 
that can help determine the cause of cosmic acceler- 
ation in the near future can be addressed most eas- 
ily. The DGP model itself is strongly constrained by 
current observations, see e.g. [15, 59-62]. 

Previous work in the literature discussed gravita- 
tional instability in this model in different approx- 
imations: [47] obtained the quasistatic linear and 
nonlinear evolution in the spherical approximation, 
while [63] developed the linear quasistatic approxi- 
mation further by going beyond spherical symmetry 
and including an explicit treatment of the bulk per- 
turbations, which justified the assumptions in [47]. 
To follow perturbations of wavelength larger than 
the Hubble radius one must go beyond the qua- 
sistatic approximation; this was first done in [64] 
using a scaling ansatz and later corroborated by [65] 
using an exact numerical solution of the linear equa- 



tions. Except for [47] all these treatments work in 
linear perturbation theory using the master variable 
formalism developed in [66, 67] to solve for the bulk 
perturbations. 

In this paper we extend previous work in several 
directions. First, we derive the linear solutions with- 
out relying on the master variable formalism, follow- 
ing the standard treatment of cosmological pertur- 
bations in GR after choosing suitable perturbation 
variables based on the properties of gauge transfor- 
mations along the extra dimension [126]. Next we 
tackle the nonlinear case, focussing on the recov- 
ery of GR by the Vainshtein mechanism using a re- 
summation technique to obtain the modified Poisson 
equation propagator. Finally, we apply these results 
to the calculation of the nonlinear power spectrum 
and the corrections induced in the bispectrum due 
to the nonlinear dynamics in the modified gravity 
sector. We then comment on the impact of these 
results on looking for modified gravity in observa- 
tions. In a companion paper [48], we develop an 
N-body code and test some of the analytic predic- 
tions developed here. For similar work on N-body 
simulations see [68]. 



II. FIVE-DIMENSIONAL METRIC 
PERTURBATIONS 

We start by reviewing metric perturbations in the 
context of braneworlds, in which the observable uni- 
verse is a 3-f 1 brane embedded in a 4-1-1 dimensional 
spacetime. Most of this section contains well known 
results (see e.g. [69, 70] for reviews), but it serves to 
fix our notation, explain our unusual gauge choice, 
and motivate our strategy to solve the bulk equa- 
tions described in section III D below. 



A. SVT Decomposition 

In this paper we are interested in the scalar per- 
turbations to the metric as seen by an observer in 
the brane. Following the standard SVT (Scalar, Vec- 
tor, Tensor) decomposition (see e.g. [71, 72]) slightly 
generalized for one extra dimension, we can write in 
general for the scalar part of the metric. 



3 



e^^{l + 2(1)) dt^ + 2 V,B dx'dt 

(1 + Sg44) dw^ + 2ViC dx'dw + 2F dtdw, 

(1) 



where xA = {t. X, y, z. w) with A = 0, . . . , 4, u> is 
the coordinate on the extra (4th spatial) dimension 
with w = denoting the unperturbed brane posi- 
tion, and a{t, w) and ^{t, w) characterize the back- 
ground solution 



a + £ \w\ a. 



l + e\w\ 



(2) 



where a{t) is the scale factor, and dots denote time 
derivatives. We assume symmetry, all metric per- 
turbations except C, F are functions of C, F in- 
stead arc odd (henceforth we work with w > 0). In 
the absence of perturbations the metric in Eq. (1) re- 
duces to the background solution found in [51], with 
e = corresponding to the self- accelerated branch 
and £ = — 1 to the normal branch. From now on 
we will concentrate on the self-accelerated branch, 
although occasionally we will discuss how some re- 
sults change in the normal branch. Note that, al- 
though not obvious from Eq. (1), the backgroimd 
bulk spacetime {w > 0) is Minkowski, by making a 
change of coordinates this can be made explicit but 
the brane will have a complicated trajectory in that 
coordinate system, see e.g. [67]. Finally, we assume 
a spatially flat universe throughout this paper. 

A few comments are in order regarding the form 
of Eq. (1). The SVT decomposition implied there is 
based on the behavior of metric perturbations with 
respect to transformations of the 3 spatial dimen- 
sions of the brane. A bulk observer, could perform 
a similar SVT decomposition but with respect to 4 
spatial dimensions, which leads to 4S, 6V and 5T 
modes that would correspond respectively to spin 
0,1,2 in a local incrtial frame. The 5T modes are 
the true physical degrees of freedom (only ones that 
propagate in the absence of sources) for a spin 2 
massless graviton in 5D, the rest are induced by 
stress energy fluctuations. Out of the total of 15 
modes, however, 2S and 3V modes can be set to zero 



by appropriate choice of gauge, leaving 10 physical 
modes (2S, 3V and 5T). However, because the back- 
ground spacetime does not have symmetry under full 
4D spatial rotations, these SVT modes do not evolve 
independently. Instead, it is easier to work with the 
3D SVT decomposition modes which do evolve inde- 
pendently thanks to the fact that w = const, hyper- 
surfaces in the background space are homogeneous 
and isotropic. The price to pay is that some V and 
T modes in 4D project into S modes in 3D. That's 
why apart from the usual 4 scalar modes encoded in 
4>, ip, A, B, wc have to deal with 3 new scalar modes 
described by C, F, dg^. 

The 3D SVT decomposition divides metric pertur- 
bations into 7S (shown in Eq. 1), 6V and 2T modes. 
Gauge transformations can be used to set 3S and 
2V modes to zero, leaving 4S, 4V and 2T physical 
modes (again, a total of 10 physical modes). For 
example, the scalar modes B and F in Eq. (1) cor- 
respond to linear combination of 4D S and V modes 
(spin and 1). The scalar mode C corresponds to 
a linear combination of 4D S,V and T modes (spin 
0,1,2). We now discuss briefly how to use suitable 
gauge transformations to eliminate 3 out of 7 scalar 
modes. 



B. Gauge Transformations and Gauge Choices 

Gauge transformations for scalar perturbations 
can be written as 



5x' 



5x^ = {6x'^,Vi5xs,5x'^). (3) 



These three functions can be chosen to eliminate 

some of the scalar potentials in Eq. (1), keeping 
in mind that there is an additional degree of free- 
dom we have not yet discussed, corresponding to the 
brane position [73] which will not stay fixed at w = 
in the presence of stress energy fluctuations. To see 
which scalar metric perturbations can be set to zero 
we look at the transformation properties of metric 
perturbations under a scalar gauge transformation, 



ip ^ tp — a Sx^ — a' Sx'^, 



(4) 



(5) 
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A + 2e^°' 5xs (6) 

B^B + e'^^5x^ -e^^dxs (7) 

C^C-5x'^- e^" {5xs)' (8) 

F^F + e^^ {5x^)' - 5x^ (9) 

Sga ^ Sga + 2 (<5a;^)' (10) 

where dots and primes denote time and w deriva- 
tives, respectively. 

Prom these equations one can deduce a few useful 

gauges. The simplest to sec, is a generalization of 
the longitudinal gauge, in this case one sets 

A = B = C = 0. (11) 

This can always be obtained, since these conditions 
uniquely fix 5xs,5x'''' ^Sx'^ , respectively. An impor- 
tant property of this gauge is that the brane is bent, 
i.e. the position of the brane is not at w = in the 
presence of fluctuations, but rather at some function 

W = Wb{x^). 

The gauge that we use in this paper can be easily 
obtained from the longitudinal gauge, Eq. (11). We 
make a gauge transformation w w — Sx"^ that only 
involves the fourth spatial coordinate (thus keeping 
A = B = 0), such that (fe^)' is fixed by the condi- 
tion Sg44 = and the w-independent part of 6x^ is 
fixed by the condition that the brane is straight in 
this gauge, i.e. wt — 0. This can always be done be- 
cause it is possible to set to zero the brane position 
with a w-independent shift of coordinates. There- 
fore, our gauge reads 

A = B = Sga = = (12) 

Note that in going from the longitudinal gauge to 
our gauge we generate a nonzero C which is given 

by, 

1 r 

C=-y^ 5gii{xi',\w\)dw-\w,,{xi')\\w\', (13) 



where quantities with hats are in the longitudinal 
gauge, and the form of the second term respects 
that C must me odd in w while bringing back the 
brane to the unperturbed position. Close to the 
brane, say w — Q'^ , the first term in Eq. (13) does 
not contribute and C is directly given by the brane 
bending W{,(a;''). Thus we may call C the brane- 
bending m,ode. In geometric terms, V^C determines 
the brane extrinsic curvature at subhorizon scales. 

Since our gauge can be imposed at any w, all 
derivatives of A, B, (^944 with respect to w will van- 
ish as well. In the following we present the equations 
of motion always in this gauge. Note also that the 
induced metric at the brane is given in terms of the 
five-dimensional metric by g^,y{t,x^) = [g'll'J]w=o- 

The most popular choice in the literature is the 
so-called Gaussian Normal (GN) gauge, where one 
fixes 6x^ as in Eq. (12), but instead fixes {Sxs)' from 
the condition C = and (Sx'^)' from the condition 
F = 0. To fully specify the gauge one fixes the w- 
independent part of Sxg by setting A = at the 
brane and that of Sx'^ by setting B = at the brane. 
Since these two choices can only be taken at w = 0, 
the tw-derivatives of A and B in this gauge will not 
be zero at the brane (and elsewhere), and one must 
keep track of them. In this sense the GN gauge seems 
a bit cumbersome, although it is perfectly fine as a 
gauge choice [127]. 

C. ADM form 

To understand the role of the brane-bending more 
C in geometric terms, see e.g. [53, 54], it is useful to 
cast the metric in ADM form [74], 

ds'^ = g^u {dx" + N" dw){dx'' + N" dw) - {N dwf, 

(14) 

where g^v is the 4D metric of a hypersurface at fixed 
w and from Eq. (1) the lapse N and shift Ni^ func- 
tions are given by, 

7V = l + ^, N^ = {F,VC), (15) 

where we have used linear perturbation theory for 
N . It is also useful to write the extrinsic curvature 
tensor of a hypersurface at fixed w, 
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(16) 



where the covariant derivatives are with respect to 
the 4D metric g^,^. From this we have, for the back- 
ground metric (bars denote background values) 



^00 



1 



{e'n'Sij, (17) 



where primes denote derivatives in the extra spatial 
dimension. Therefore, the extrinsic curvature scalar 
reads 



K=-{3a' + (3'), (18) 
whereas for perturbations we have to linear order 



SK 



-2/3 



F{3a-I3) + F +{3^'-(f>')-e-^"V^C. 

(19) 

A similar (though longer) expression holds for the 
full extrinsic curvature tensor. This will imply (af- 
ter we solve the bulk equations and develop the qua- 
sistatic solution) that the extrinsic curvature will 
be dominated by second derivatives of the brane- 
bcnding mode at subhorizon scales, and by normal 
derivatives at larger (but still smaller than Hub- 
ble) scales. At scales larger than the Hubble radius 
the contributions from F become important, but we 
won't deal with such wavelengths in this paper, thus 
F will play essentially no role. Since extrinsic curva- 
ture is the new degree of freedom that makes gravity 
different from GR, this shows which metric pertur- 
bations are going be relevant for structure formation. 



where GJ^g is the 5D Einstein tensor, and Gab, Tab 
the standard 3-1-1 Einstein and stress-energy ten- 
sors, which vanish when any index is larger than 
3; this condition for Tab enforces there is no flux 
of cncrgy-momcntiim to the extra dimension (recall 
the position of the brane is fixed at w = even in the 
presence of fluctuations). In Eq. (20), g^^^ denotes 
the absolute value of the 5D metric determinant, and 
[Su] = Sr){'w)^/g^. The parameter Vc denotes the 
cross-over scale from 4D to 5D behavior, indeed one 
can rewrite by dimensional analysis Eq. (20) for a 
point source at the brane as 



1 1 



rr. 



(21) 



where cj) denotes a generic component of the metric. 
Therefore, for r <C rg one has the standard 4D be- 
havior <^ ~ r~^, whereas for r ^ Vc the 5D behavior 
(f) ^ r^^ is obtained. More precisely, small pertur- 
bations about Minkowski ^j^^ = ry^i^ -|- hixi/ 

(V < 1) 

the metric at the brane can be readily obtained in 
momentum space [76] 



G{p) (V-3»?^.r) 



= G{p) (f^, - ^7?^, f ) + ^G{p) rj^^ f, 

(22) 

where a tilde denotes objects in Fourier space (de- 
pending on four- momentum p'^),T = T^, and 



G{P) 



inG 



(23) 



III. FIELD EQUATIONS AND LINEAR 
THEORY 



A. Basics 

The field equations corresponding to brane- 
induced gravity in five dimensions are [75] , 



^^'^ Gab + Gab [fo] = SttG Tab [fo], (20) 



2r, 



from which it follows Eq. (21) for a static point 
source. We will come back to Eq. (23) when we 
discuss linear cosmological perturbations as the un- 
usual p^^ dependence will lead to interesting scale 
dependence of the growth factor at large scales. In 
the second line of Eq. (22), we have decomposed 
the tensor structure of the full theory (previous line) 
in terms of the standard massless spin-2 part (first 
term) and the helicity-zero part (second term) that 
effectively acts as a scalar. The fact that the full 
theory in the massless limit pre ^ 1 does not re- 
cover the massless case (GR) is the well-known van 
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Dam-Veltman-Zakharov (vDVZ) discontinuity [77- 
79]. As conjectured by Vainshtein in massive grav- 
ity [43] , the discontinuity should disappear in the full 
nonlinear theory (in which the helicity-zero mode 
gets suppressed) , otherwise one could rule out a mas- 
sive graviton (m ~ here) from local observations 
no matter how tiny its mass. Understanding this 
Vainshtein mechanism for realistic cosmological per- 
turbations is one of the main goals of this paper. We 
shall see that the behavior of the graviton propaga- 
tor jTist discussed has a close analogy to the modified 
Poisson propagator for cosmological perturbations, 
as expected. 

It will be useful for future reference to write the 
gravitational part of the action for Eq. (20) in ADM 
form, 



1 



IGttG 



1 



+ --N{R + K^ -K^^K'^") , (24) 



(28) 



where = dx^/dr = {a,v^) with v the standard 
velocity fluctuations about the Hubble flow, and t 
is the conformal time. The only non- vanishing com- 
ponent of the stress tensor fluctuations are (in linear 
perturbation theory), 



6T°=pS, 
ST^ = —apvi 



(29) 



(30) 



From this, conservation of energy and momentum 
gives respectively (always to first order). 



a 



(31) 



where R is the Ricci scalar for the 4D metric. The 
first term is the usual Einstein-Hilbert term, whereas 
the second contribution proportional to the lapse 
fimction N is the 5D Einstcin-Hilbcrt term decom- 
posed in ADM form into intrinsic (i?) and extrinsic 
{K) curvature of hypersurfaces at constant w. For 
the Ricci scalar of the background we have, 



R = -Qe-'^^{2a^ - d/3 -h a) 
whereas for perturbations. 



(25) 



6R = 2e-2"V^(^-2^/>)-2i?(/- 

+Qe-'^'^{a4) + 4aip - ^ip + i)) (26) 



a 



(32) 



The first is the standard continuity equation with 
the extra term representing the change in volume 
as the universe expands due to the fact that a(l — 
tp) is the perturbed scale factor. The second is the 
geodesic equation for the fluid, in the approximation 
of weak fields and small velocities. These equations 
are independent of the theory of gravity (as long as 
it is a metric theory and stress energy is conserved), 
and must be supplemented by the field equations of 
the gravity theory under consideration, which relates 
stress energy fluctuations to metric fluctuations. 



B. Background Evolution 



The stress-energy tensor we assume throughout 
this paper is that corresponding to a dark matter 
fluid. 



(27) 



Let's briefly recall the evolution of the background 
in the self-accelcratcd branch, this will also be useful 
to set the time variables we are going to use when 
discussing perturbations. From the 00 fleld equation 
we obtain the modifled Friedmann equation, 



where p is the mean density, 6 describes density fluc- 
tuations and the four-velocity. 



H^ = l.l.G-p, 



(33) 
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Since p > 0, the theory makes sense only when 
rcH > 1, thus probing distances larger than the 
crossover scale rc necessarily implies considering per- 
turbations of wavelength larger than the Hubble ra- 
dius. As time goes on, one immediately sees from 
Eq. (33) that since p ^ 0, one enters a de Sit- 
ter phase where the universe expands exponentially 
with H = l/vc- The acceleration follows from the ii 
field equation, which gives 



TrH-l 



fj — — TJ 

2 reF-1/2' 



(34) 



which goes from —1/2 in the early universe to 1 
as the universe gets into the purely de Sitter phase. 
We can then write simply. 



e"" = a{l + \w\H), el^ = l + \w\Hq. (39) 
Finally, note that before the universe enters into 
the acceleration phase (77 > 2), g < and from 
Eqs. (1) and (39) it follows that there is a Rindler 
horizon in the bulk where (3 ^ —00 located at 



which implies that the acceleration of the Universe 
is given by 



Wh 



-iqH)- 



(40) 



a 2r] 



1' 



where we defined, 



(35) 



(36) 



thus the Universe starts to accelerate when H drops 
below the critical value r/crit = 2. In the early uni- 
verse T] ^ 1, in fact T] can be thought of as a time 
variable which decreases with cosmic time (77 > 1, 
achieving = 1 in the infinite future when the dark 
matter density asymptotes to zero). In terms of this 
time variable, the scale factor reads. 



viv - 1) -I 



1/3 



(37) 



and rjQ = H^rc can also be written in terms of the 
2 = matter density, r/o = (1 — fij^)"^. One can go 
back and forth between o and r} using Eq. (37) and 
its inverse relation, 277 = 1 + -^/l + 4:i]a{ria — l)o~^. 
We can define the acceleration parameter. 



aa 



27] -r 



(38) 



which becomes 2H^^ in the early universe, and goes 
to infinity when the universe approaches the accel- 
erated phase, and disappears afterwards. 

In what follows we will use rj and q to character- 
ize the time dependences of the scale factor in this 
model. 



Linear Perturbations: at the Brane 



We now consider linear metric perturbations at 
the brane. Since the dark matter has negligible 
anisotropic stress, STj = 0, and thus from the in- 
equations {i ^ j) a.t w = one obtains the follow- 
ing simple boundary condition for the brane-bending 
mode C, 



Co = rc (V'o - <l^o), 



(41) 



where a subscript denotes quantities evaluated at 

the brane. The OO-cqiiation at the brane gives the 
modified Poisson equation in linear theory. 



(42) 
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At this point it is worth making a few observations. 
First, note that we can obtain the GR hmit by tak- 
ing ijjQ = (po in Eq. (41) which sets to zero Co 
in Eq. (42), and then ignoring the contribution in 
parenthesis (effectively setting Vc to infinity) gives 
the Poisson equation in GR. The terms in braces on 
the right hand side of Eq. (42) correspond to the 
usual contributions in GR from retardation effects 
that dominate at super horizon scales, while at sub- 
horizon scales they are negligible compared to V^tpo- 
The reason that normal derivatives only appear at 
most to first order in Eq. (42) is that at least an- 
other first normal derivative of the background met- 
ric (that has a jump at w = 0) is nec^ded to got 
a nonzero contribution localized in an infinitesimal 
brane. 

We are interested in finding the growth of struc- 
ture at subhorizon scales where time derivatives are 

small compared to spatial gradients, and thus one 
can assume a quasistatic evolution. In this case, 
Eq. (42) becomes, 



^V'(.^o+V'o) =47rGa2p5, (43) 

where we have used the constraint in Eq. (41) to re- 
place the brane-bending mode Co in terms of poten- 
tials ijjQ , (j)Q , and have neglected the terms in paren- 
thesis and braces as they involve spatial derivatives 
of order less than second supressed by the Hubble 



radius or 



Then, all that remains to find 



is the relation between 00 and tpQ. This can be ob- 
tained directly from the 44 bulk equation close to the 
brane {w = 0"*"), which in the same approximation 
reads. 



V2(Vo - .^o) = 



1 



viq + 2) 



V2(2Vo - (44) 



and thus the system is closed. This is essentially 
the argument given in [47] to derive the growth fac- 
tor at subhorizon scales. Note that Eq. (44) relates 
the extrinsic curvature, Eq. (19), to the Ricci scalar, 
Eq. (26), at subhorizon scales; we shall come back 
to this below. At high redshift, r] = Hvc ^ 1 while 
q -1/2, and Eq. (44) sets Vo = <f>o and GR is 
recovered in Eq. (43). 

The modified Poisson equation that results from 
Eqs.(43)-(44) is. 



VVo = 47rGeffa2p(5, G^s = G 



2r,iq + 2) 



2?7(9 + 2)-3' 
(45) 

thus modification of gravity in this approximation 
can be understood as a time-dependent (but scale- 
independent) modification of the gravitational con- 
stant. Note that at late times, in the de Sitter limit, 
77, g ^ 1 and Gcff (2/3)G. Thus gravity becomes 
weaker leading to a slower growth of structure com- 
pared to a dark energy model with the same expan- 
sion history. 

These arguments, while certainly reasonable, are 
not rigorous enough to understand the limit of ap- 
plicability of the approximations made. For exam- 
ple, that Gcff is only a function of time results from 
keeping only the highest (second) spatial derivatives. 
We would like to understand the limitations of these 
results at both large and small scales. The limita- 
tions at large scales come from neglecting the nor- 
mal derivatives to the brane, ipQ, for which we need 
to solve the equations in the bulk [63]. The limita- 
tions at small scales come from the fact that brane- 
induced gravity becomes GR at small scales through 
the Vainshtcin mechanism, where nonlinear effects 
suppress the brane-bending mode Co. For this we 
need to understand the nonlinearities responsible for 
this mechanism and study at what scale they do op- 
erate. We first consider the linear bulk equations, 
postponing until section IV the discussion of nonlin- 
ear effects. 



D. Linear Perturbations: in the Bulk 

Before we discuss the solution of the bulk equa- 
tions, it will be useful to split perturbations into 
two classes according to their behavior under gauge 
transformations that shift the brane position, w — > 
w + 5x'^. We define two new metric variables, 

$ = + ^ + (/J' _ a') C, f = F-C, (46) 

which, as can be seen from Eqs. (4-10), are in- 
variant under shifts of the brane position. We also 
define a new metric variable, 



* = 2V' - 



(47) 
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which gives the leading contribution to the per- 
turbed Ricci scalar, Eq. (26), at subhorizon scales. 
We keep the brane bending mode as our fourth met- 
ric variable, thus we are looking for solutions in 
the bulk for $, /, vj/, (7, where the first two variables 
are invariant under gauge transformations involving 
only w, and the other two are not. As we shall see, 
this split between <I>, / and C will play a key role 
in solving the bulk equations. 

The equations that determine the shift vector 
in the bulk are very simple. The ij field equation 
reads, 



in some cases rather complicated but wc are in- 
terested in sub-horizon perturbations where time 
derivatives are suppressed compared to spatial gra- 
dients. Therefore, in the following we make the ap- 
proximation that wc consider modes inside the Hub- 
ble radius, i.e. {k/aH)"^ 1, where k is the comov- 
ing wavenumber. This is a very good approximation 
since, e.g. at ^ = 

_|-=300 f .y (53) 

aH Vo.l/iMpc"^/ 

In this approximation, Eq. (50) becomes 



3C" + 2{2a' + p')C = 2^- 



(48) 



which determines C once $ and \E' are known. After 
using this result, the Oj field equations give 



/'+(a'+/3')/ = 2($ + a$). 



(49) 



This determines / from the evolution of An equa- 
tion for $ follows from subtracting one third of the 
ii field equations from the 00 equations. 



□5 $ = 4de-^^$ + 4(a' - /3')*' + /, C}, (50) 

where the 5D d'Alembert operator in the back- 
ground metric reads 



□5 = □ - (3a' + /?') dy, - dl 



(51) 



and the 4D one, 

□ = e-2^ d'l + (3d - /?) e-2^ dt - e-2° V^ (52) 

and the term in braces in Eq. (50) denotes terms 
that involve the four metric variables (but none of 
their derivatives). All the terms that do not involve 
$ have the property that their coefficients vanish in 
the de Sitter limit (g = 1), as can be seen explicitly 
for the term, i.e. a' = = -ff/(l + wH) in 
the de Sitter phase. In this limit, {} = Id? e~'^^ 
and Eq. (50) becomes a closed equation for <I>. This 
shows that during the accelerated phase we are very 
close to having decoupled the bulk equations, even 
beyond the quasistatic approximation. 

So far we have written exact bulk equations in 
linear theory. The remaining bulk equations are 



$" + (3a' + /?')$' + c"2«v2$ 



0, (54) 



where we have also neglected 4*', / and C terms 
(all of which disappear as the universe approaches 
the de Sitter phase). This will be further justi- 
fied shortly. Equation (54) can be readily solved 
in Fourier space as 



$ = $0 (1 + s)-'=/«^ s = wH 



(55) 



where we have dropped the solution that diverges at 

infinity. Strictly speaking, before the universe accel- 
erates one must impose vanishing boundary condi- 
tions at the Rindler horizon given by Eq. (40), but 
the quasistatic approximation breaks down at such 
large distances away from the brane, and in any case 
Eq. (55) decays so fast that any mixing with the sec- 
ond solution is going to be negligibly small. Equa- 
tion (55) is analogous to the quasistatic bulk solution 
obtained in [63] for the master variable (however, it 
is not the same as the master variable is a gauge in- 
variant quantity, and $ is only invariant under gauge 
transformations involving only the extra-dimension 
coordinate). 

We can check the validity of the quasistatic solu- 
tion by comparing the time derivatives to the other 
terms in the Eq. (50). We have 



1 + s 



q ln(l + s) 



(56) 



and compared to spatial derivatives in the left hand 
side of Eq. (50) we see that time derivatives become 
important only far from the brane, i.e. when s ~ 1 
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or ui ~ H~^; at this point the qiiasistatic solution 
in Eq. (55) has decayed exponentially to zero for the 
wavenumbers of interest, see Eq. (53). 

Having a solution for $ we can now solve for / 
from Eq. (49). Here one must be careful to include 
time derivatives that act as a source of / in Eq. (49), 
leading to 



f' 



{l-q)s+q{l+s) ln(l+s) 



2e'«V-2$', 
(57) 

where we have set the homogeneous solution /h oc 
{l+s)~^{l+qs)~^ to zero since it violates the i4 bulk 
equations, as we shall see in Eq. (65) below. Note 
that close to the brane a sublcading term of order 
$ {aH/k) becomes dominant since Eq. (57) vanishes 
as s ^ 0. At the brane (s = 0+) we then have 



In what follows we neglect this sublcading contri- 
bution, and use Eq. (57) everywhere in the bulk. 
Then / is at most of order Note that neglecting 
/ altogether results in a shift vector that is a pure 
four-gradient, i.e. A^^ = V^C 

It remains to find the dynamics of ^ and of the 
brane-bending mode C which determine, respec- 
tively, the Ricci scalar and the extrinsic curvature in 
constant w hypersurfaces at subhorizon scales. As it 
is clear from the ADM formalism, Eqs. (15) and (24), 
the 44 field equation relates the Ricci scalar to the 
extrinsic curvature (the Gauss-Codazzi equation). 



R = K^^K'"' - K\ 



(59) 



/o « -2$. 



which to linear order gives, in the subhorizon ap- 
(58) proximation (using the same notation as in Eq. 50) 



{2a' + f3')DC - + (/?' - 3q!')*' = e'^'^ l> - f3'^' + {/, C, C, 6, 



(60) 



Again, working in the quasistatic approximation, 
and neglecting ^I*' compared to V^^, which will be 
justified shortly, we obtain a relationship between * 
and C that can be written, 

* = (2a'+/J')C+e^'"-''V-2|.-/J'e2«V-2$'. (61) 

In principle one could add to this the solution of 
□^1 = (/?' - SaO^-', which is * = (1 + s)-^'/3(l + 

2qs/{3 — q))^^ with k = k/aH, but such term 
would violate the i4 field equation below. This is 
important because otherwise the solution for the lin- 
ear growth factor would be different at subhorizon 
scales. Replacing Eq. (61) in Eq. (48) and ignoring 
subleading terms gives a very simple equation for C 
in terms of $ alone, 



3C" = 



which leads to 



C = Co + ^(^)[(1 + «)-'=/»« -1 



(62) 



(63) 



and this in turn gives a solution for ^ from Eq. (61). 

We are now in a position to check the validity of 

the approximations made. First, to go from Eq. (50) 
to Eq. (54) we have dropped the terms involving /, 
^, ^' and C, apart from assuming the quasistatic 
approximation for $. Since / is at most of order $, 
Eq. (57), it is justified to drop it compared to V^$. 
The remaining terms involving 5", VE"' and C are po- 
tentially dangerous, as C does not decay into the 
bulk and ^ only decays slowly, thus they can become 
larger than as it becomes strongly damped in 
the bulk. However, using Eqs. (61)-(63) it follows 
that all these terms in Eq. (50) cancel leaving only 
contributions proportional to $. This is as expected 
from the properties of fields under gauge transforma- 
tions that shift the brane position. Finally, we see 
from Eqs. (61)-(63) that it is also justified to drop 
compared to V^^* in the quasistatic regime. 

Other checks can be done from the remaining bulk 
equations. The i4 field equation reads. 
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f + {a-p)f = --e^'^ 3*' + (a'-/3')(2$-* 



+ (8a'^-a'/3' + 2/3'2)C 



(64) 



which in view of Eqs. (61-63) can be written as 



(65) 



which is solved by Eq. (57), and does not support 
a slow decay into the bulk, i.e. the homogenous so- 
lution of Eq. (49), /h oc (1 -h s)"Hl + '7s)"\ does 
not solve Eq. (65) with zero right hand side, and 
thus must be discarded. Finally, the 04 bulk equa- 
tion can be shown to hold by using Eq. (57) for / as 
well. 

These results are summarized in Fig. 1, where 
we show the four bulk solutions as a function of 
the coordinate into the bulk in units of the Hub- 
ble radius, s = wH, for two different wavectors 
k = 0.005 ft. Mpc"^ (solid) and k = 0.02ftMpc"^ 
(dashed). Note that fields invariant under shifts of 
the brane position ($, /) are concentrated at the 
brane and leak only partially into the bulk, whereas 
those who are not (C, ^) approach a constant far 
from the brane. However, one should keep in mind 
that this neat separation of bulk behavior according 
to symmetry breaks down at small scales (high-fc) 
when nonlinearities are included, as we shall discuss 
in section IV. 




1 1 ■ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ■ 1 1 ■ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 i 1 1 1 1 1 1 1 1 1 1 1 
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
s=wH 



FIG. 1: Bulk solutions, for k = 0.005 /iMpc~^ (solid) 
and k = 0.02/iMpc"^ (dashed), as a function oi s = wH 
for redshift z = 0, when the matter density is given by 
0,%^ = 0.2. The brane is located at s = 0. All metric 
variables are normalized by the value of $ at the brane, 
$0- Fields invariant under shifts of the brane position 
($, /) decay strongly into the bulk, whereas C, ^' do not. 



E. The Linear Quasistatic Solution 

Now we can go back to the brane equations 
and finalize the linear quasistatic solution. From 
Eqs. (46)-(47) we can map back to the (j>.ip,F vari- 
abls relevant to the Poisson equation, Eq. (42), 



2(a'-/?')C + 2$-* 



(a' - I3')C + ^ + ^ 



(66) 
(67) 



from which we can obtain the normal derivatives at 
the brane. 



00 ~ 2Vjo 



2k 
( 

3 a 



(68) 



plus subleading (order i?$o) contributions. We thus 
see that normal derivatives are suppressed compared 
to the leading contributions in the Poisson equation 
Eq. (42) by one inverse power of krc/a. These terms 
are genuinely of quasistatic nature (i.e. they give 
rise to the 5D decay of the gravitational potential 
for static sources), and thus we must keep them as 
part of the quasistatic solution. They are dominant 
compared to retardation effects standard from GR 
(which are also modified here at late times by 5D dy- 
namics) , which are suppressed by two inverse powers 
of k/aH. Clearly, krc/a compares a physical scale 
a/k against the spatial scale rc, whereas k/aH com- 
pares the light crossing time across a perturbation 
of size a/k to the Hubble time a measure of 

retardation effects. 

An important feature of brane-induced gravity 
is that from the point of view of standard 4D 
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physics, the bulk dynamics introduces tlirough nor- 
mal derivatives such as those in Eq. (42), unusual 
non-local operators. From Eq. (68) we have that 
normal derivatives of $o can be written in configu- 
ration space as 



third goes as 0/(rrc), then we can symbolically write 
the Poisson equation as (compare to Eq. 21), 




$^ = --yZ^$o. (69) 
a 

This is a consequence of the 5D nature of the gravity 
theory. 

We can then summarize the linear quasistatic 
equations as follows. To make equations compact, 
we will first make use of different (though linearly 
related) metric variables. We later write down every- 
thing in terms of the Newtonian potential (?io- Using 
Eq. (61) at w = 0+ we have, 

H{q + 2)V^Co + Hqa^/^^o = ^^^o, (70) 

which is a more accurate version of Eq. (44). To 

write the Poisson equation, Eq. (42), in the qua- 
sistatic approximation, we can drop the terms in 
braces in Eq. (42) (as usual in GR) since they are 
down from the leading terms by the factor (aiJ/fc)^, 
and also we can drop Fq as from the bulk solution 
— Co, and thus this term is subdominant by the 
same factor to the term with V^Cq. Then we have, 

VVo-;^V2Co-:^\/-V2$o = 4nGa''pS. (71) 

Note, as argued above, that the nonlocal terms 
in these equations belong to the quasistatic approxi- 
mation. For a point source of mass m, (5 ^ l/^'', 
and the right hand side of the Poisson equation 
scales as Vg/r^, where oc Gm is the gravitational 
(Schwarzschild) radius of the source. The first two 
terms in Eq. (71) scale as where is the or- 

der of magnitude of metric perturbations, while the 



and thus the normal derivative term governs the 
transition of the static on-brane potential from the 
usual 4D behavior cf) ~ I'g/f for r < Tc to the 5D 
behavior (j) ^ [Tg^c)/''''^ for r > Tc- Eq. (71) thus 
reflects precisely the behavior of the propagator in 
Eq. (23) for static sources, with the nonlocality in 
the Poisson equation a consequence of the be- 
havior of the on-brane propagator. This is a feature 
of extra-dimensional theories of gravity that cannot 
take place in a local model of dark energy (no matter 
how complicated), and it is interesting to note that 
it is expected to happen in theories with more than 
one extra dimension as well, as long as there is a 5D 
regime before hitting the truly large distance regime 
(as in "cascading" models, see [36-38]). 

As we shall see, these non-local corrections be- 
come interesting well before scales of order Tc and 
improve the validity of the quasistatic approxima- 
tion at large scales compared to the quasistatic so- 
lution of [47, 63] which sets the normal derivatives 
to zero. Therefore, it is expected that at least some 
of the deviations reported at large scales in numeri- 
cal solutions that include normal derivatives [64, 65] 
can be explained by this simple result. 

Note that in Eqs. (70)-(71), the brane metric vari- 
ables are all related through Eqs. (41), (46) and (47). 
Thus there is only one independent brane potential 
that obeys a Poisson- like equation. We want to com- 
pute the growth of density perturbations, and since 
the Newtonian potential is responsible for accel- 
erations, Eq. (32), we write everything in terms of 
the Newtonian potential. Thus we have the effective 
Poisson equation in Fourier space 



4nGa'^p6, or - k'^(j)o = 4TrGes{k,ri)a^p5, (73) 

where we have defined an effective gravitational constant GcS that depends on scale and time [128]. For 
completeness, we note that one can also derive the following Poisson-like equation for the brane-bending 
mode. 



_ 2 , f 27?(2 + g)-3 , /3a v r/(l + g)-1 ^2^ _ 
l2r?(2 + g)-4 + UreA,,(2 + g)-2J /" 
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-k' ^ {[2,7(2 + g) - 3] - (^)(2W - 1)^ + q) - 2]} ^ SnGa^pS. (74) 



oO.B 



0.6 























:/ z=o 





0.001 0.01 0.1 1 



k [h/Mpc] 

FIG. 2: The effective gravitational constant Gcft{k,z) 
that governs the hnearized Poisson equation for the New- 
tonian potential, Eq. (73), as a function of k for different 
redshifts z = 0, 0.5, 1, 2, 3, 10 (from bottom to top). Non- 
linear effects modify significantly this behavior at small 
scales, see Fig. 6 below. 



Figure 2 shows Ges{k,z) as a function of k for 
different redshifts. We see that at all times GcS is 
less than G, more so at late times (when the modi- 
fications away from GR become more important, as 
the universe accelerates), and at large scales where 
the transition to the 5D regime of even weaker grav- 
ity takes place, inducing scale dependence. For the 
purpose of this plot we have assumed fi^j — 0.2 for 
which 770 = 1-25 and the universe starts accelerating 
at z = 0.86. 

The two features seen in Fig. 2 are easy to under- 
stand from the basic physics of branc-induced grav- 
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ity. The scale-indepedent suppression of GcS arises 
from the first term within braces in Eq. (73), which 
is clearly always larger than unity, while the scale de- 
pendence arises through the second term that is al- 
ways positive, enhancing the previous effect at large 
scales. 

The physical origin of these terms is very simple. 
They arise because the gravitational effect of sources 
at the brane is not only to create 4D intrinsic cur- 
vature (encoded by the 4D Ricci scalar), but also 
to create extrinsic curvature, which at small scales 
is dominated by the brane-bending mode Cq (first 
term within braces in Eq. 73) and at larger scales 
by normal derivatives (second contribution within 
braces in Eq. 73), see Eq. (19). 

The precise fraction that goes into the usual 4D 
Ricci scalar (available to accelerate dark matter par- 
ticles and thus grow structure) versus extrinsic cur- 
vature is controled by the Gauss-Codazzi equation, 
Eq. (59), that relates intrinsic to extrinsic curva- 
ture. Therefore, as long as the extrinsic curvature 
remains important, density perturbations will have 
suppressed clustering because part of their ability to 
induce intrinsic curvature (the only type available in 
GR) is being used to create extrinsic curvature of the 
brane. We will see in the next section that nonlin- 
ear effects suppress the extrinsic curvature at small 
scales leading to the recovery of GR. 

From the modified Poisson equation, Eq. (73), and 
the conservation equations, Eqs. (31-32), at subhori- 
zon scales we have the linear density perturbation 
growth equation, 

^+n^^ V^^o = 47rGeff (fc, 7^)a^pd (75) 
dr^ or 

where r is conformal time {adr — dt) and Ti = aH is 
the conformal expansion rate. Changing from con- 
formal time to our time variable, drj = T-l{q— l)r] dr, 
we can find the growth factor as a function of time 
and scale. The results are shown in Fig. 3. The over- 
all behavior of the growth factor is expected given 
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FIG. 3: The growth factor in the quasistatic approxima- 
tion as a function of k for different redshifts z — 0, 0.5, 1 
(from top to bottom). 



the properties of GcS shown in Fig. 2, the scale de- 
pendence is much smaUer due to the fact that the 
growth results from the integrated history of the 
strength of gravity and for most of the evolution the 
scale dependence in GeS is small. Our results are in 
good agreement with those in [64, 65] for the scale 
dependence at large scales, suggesting that most of 
the effect is due to the normal derivatives in the qua- 
sistatic approximation. 



IV. BEYOND LINEAR THEORY 

A. Nonlinearities 

We now discuss how the linear results derived in 
the previous sections change when nonlinear effects 

I 



are included. The first step is to identify the sources 
of new nonlinearities not present in GR. Looking at 
the action in ADM form, Eq. (24), the new con- 
tributions are in the second term proportional to 
the lapse function N. We see right away that new 
relevant nonlinearities appear through the 44 field 
equation, Eq. (59), where the extrinsic curvature is 
SK cx V^C, see Eq. (19). Since V^Cq oc (5 in lin- 
ear theory, clearly the quadratic terms in Eq. (59) 
will become important as 6 becomes of order unity. 
Neglecting subleading terms the nonlinear version of 
Eq. (60) becomes 



-2a 

(2a' -f- /?') V^C + ^-[(V^C)^ - (V,yC)2] = V^*, 

(76) 

which agrees with the equation written down in [80] 
for the brane-bending mode; see Eq. (87) below for 
a version of this equation at the brane including 
normal derivatives. Therefore the dynamics of the 
brane-bending mode is nonlinear [44, 47, 53, 54, 80- 
82]; this is crucial as it determines the relationship 
between the Newtonian potential 4>o and the spatial 
curvature potential -001 through Eq. (41), which it 
is easy to check does not receive significant correc- 
tions even in the nonlinear regime. These two facts 
together determine how the theory reduces to GR at 
small scales, we'll explain this in more detail in the 
next sections. 

There are other, similar, nonlinear corrections 
that affect the bulk solutions derived in the pre- 
vious section. The main reason is that while <I> 
obeys an approximately decoupled equation from C 
at the linear level due to its properties under gauge 
transformations involving the extra dimension, when 
nonlinear effects are included this result no longer 
holds, as fields of different symmetry are coupled. 
Thus there is a coupling between fields being ap- 
proximately constant in the bulk (^E", C) and those 
that decay strongly ($,/). For example, Eq. (48) 
becomes (fc 7^ i, j) 



V.j[3C' + 2(2a' +P')C + ^ + Se^^" [VyC VfcfcC - V,;feC V^fcC] = 0, 



(77) 
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therefore there are corrections of order {k/aH)^HC 
to Eq. (48) and Eq. (62), with the consequence that 
the details of how C behaves in the bulk will be 
modified somewhat. In second-order perturbation 
theory, is easy to see that a new contribution devel- 
ops that makes C decay slowly into the bulk as e~", 
this is expected to kick in at high k. 

Similarly, the bulk equation for Eq. (54), de- 
velops in second order a source term going as ^ 
g-2a ^^2(^-^2^ which again becomes important at 
high-fc, meaning that <1> will not decay as strongly 
into the bulk as given by linear theory, Eq. (55). 
This means that the normal derivatives at high-fc 
will not be as large as given by Eq. (68), which in 
turn implies that the nonlocal terms in the equa- 
tions of motion will go away at small scales faster 
than expected based on linear theory. Summariz- 
ing, these nonlinear effects will make small nonlo- 
cal terms even more subdominant in the nonlinear 
regime. Therefore, we conclude it is safe to neglect 
these nonlinearities in the bulk equations. 

Finally, there are other possible nonlinear terms 
that involve the interaction of C with the other met- 
ric perturbations. Since it is precisely the brane- 
bending mode C that is a ghost (at least in the de 
Sitter limit), one concern is that interaction of C 
with normal fields can develop a catastrophic insta- 
bility of cosmological solutions (see e.g. [83, 84] for 
toy models and discussion of quantum effects). 

The most dangerous of such interaction terms 
that correct e.g. Eq. (60) have the form $'V^C to 
quadratic order and $$'V^C to cubic order. Given 
that does not grow with k at high-fc due to the 
nonlinear effects discussed above and that its power 
spectrum decreases rapidly with k (P$ oc k~^), these 
terms do not give significant corrections on any cos- 
mological scales, since $ <C 1. In regions where 
the weak field approximation breaks down ($ ~ 1), 
e.g. near a black hole, the situation can of course 
be very different. Although in such high density re- 
gions, one may have difficulties exciting C due to the 
Vainshtein suppression, and since C does not obey a 
superposition principle, one may not add "free" so- 
lutions to cosmological ones, thus any analysis must 
be done self-consistently. See Fig. 7 in paper II [48] 
for a self-consistent solution showing how a cosmo- 
logical background is enough to suppress C at very 
large scales compared to astrophysical objects. 



B. Vainshtein Mechanism for Cosmological 
Perturbations 

Having identified the main new sources of nonlin- 
earities, we are ready to discuss their effect. Before 
wc do that in detail, it is useful to give a simple pic- 
ture in geometrical terms of the role of nonlineari- 
ties, which lead to the recovery of GR at small scales, 
by the Vainshtein mechanism [43]. For detailed dis- 
cussion of this in the context of the Schwarzschild 
solution see [44, 46, 81, 85-88], here we are rather 
interested in how the mechanism works in a cosmo- 
logical setting, for previous discussion along these 
hues see [47, 80]. 

We start from the field equations at the brane, 
which for 4D components can be written as 

(-^M^ - ^g^i^R) - -{K^,~g^,K) = SnGTf,,, (78) 

where the first term in parenthesis is the usual Ein- 
stein tensor characterizing intrinsic curvature, and 
the second term denotes the contribution coming 
from extrinsic curvature of the brane embedded in 
the bulk. The relationship between intrinsic and ex- 
trinsic curvatures is given by Eq. (59), and valid 
everywhere in the bulk (including arbitrarily close 
to the brane, which is what we are interested in 
this section). Roughly speaking, we can interpret 
the self-accelerated solution by saying that for this 
branch the extrinsic curvature is of the same sign as 
intrinsic curvature, leading to the Friedmann equa- 
tion Eq. (33) and a self-accelerated universe where 
intrinsic curvature (^ H^) is balanced against ex- 
trinsic curvature of the brane (~ H) when density 
of matter has dropped sufSciently. The threshold 
extrinsic curvature at which this transition happens 
is given by the inverse curvature radius r~^. There- 
fore, the universe accelerates when extrinsic curva- 
ture starts to play a role. 

We can understand the behavior of perturbations 
and the recovery of GR at small scales in similar, ge- 
ometric terms. For this purpose let's rewrite Eq. (59) 
in the following symbolic form (hereafter we also 
drop numerical factors of order unity) 

R = R + SR^ = {K + 5Kf (79) 

where quantities with bars refer to background val- 
ues. Since R = H^, for perturbations we have 
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5R^ K6K + {SKf, (80) 

with SR ~ V^*/a^ and 6K ~ V^C/o^, assuming 
we are at scales well within the Hubble radius. We 
can then write 



spherical isolated perturbations, for which 4> ~ rg/r 
where rg is the gravitational radius corresponding 
to the mass enclosed by the perturbation [47]. Since 
C ~ rc(t>, using Eq. (83) we can translate the re- 
quirement SK/K ^ 1 to find the scale away from 
the perturbation where the transition to GR takes 
place 



5R _6K /6K\-^ 



where 



R 

SR 



/ A; \2 



K \aH 



C 



(81) 



(82) 



(83) 



Then at large scales, the solution of Eq. (81) is 



5R 



5K 



(5< 1 



(84) 



and fluctuations in intrinsic and extrinsic curvature 
are comparable, thus gravity is modified by the con- 
tribution of extrinsic curvature fluctuations encoded 
in the brane bending mode C. The same is true at 
scales approaching the Hubble radius when normal 
derivatives dominate the contribution to the extrin- 
sic curvature, leading to 5D behavior. 

When 6 becomes of order unity, however, the sec- 
ond term in Eq. (81) becomes important. Therefore, 
at small scales we have instead. 



SR /SK\^ , ^ 



(85) 



and fluctuations in intrinsic curvature (~ S) are 
much larger than in extrinsic curvature (~ (5^/^). 
Thus fluctuations in extrinsic curvature of the brane 
are suppressed at small scales and GR is restored. 

Note that although fluctuations in extrinsic cur- 
vature become large in the nonlinear regime, there 
is no metric perturbation that is of order unity [86]. 
One can sec that from Eq. (83), while SK/K may 
be large, (?j4 = VjC is always small, since k/aH is 
very large in the nonlinear regime. We will explicitly 
check this with N-body simulations in paper II. 

We can put the condition for transition to GR, 
SK/K ~ 1 in a more familiar form by considering 



which gives the Vainshtein scale ~ rcH~^rg ~ 
r^rg familiar from the Schwarzschild solution [44, 46, 
81, 85- 88]. The difference here is that in the cosmo- 
logical context Vg is the gravitational radius corre- 
sponding to the mass enclosed by the perturbation, 
not the total mass including the background, whose 
gravitational effects are already taken into account 
in the evolution of the Hubble constant [47]. Indeed, 
although at small scales the relation between metric 
and stress-energy fluctuations approaches GR, per- 
turbations still evolve in a background that has a 
modified expansion history. The difference between 
total mass and perturbation mass becomes totally 
negligible for astrophysical objects within the range 
of their influence as the mass contributed by the 
background is completely negligible. For c;osmolog- 
ical perturbations, the opposite is the case for long 
wavelength perturbations and thus while counting 
total mass one would expect no corrections to GR 
below scales [45, 89] (and nonlinearities at that 
scale), properly taking into account the background 
one expects that for long-wavelength perturbations 
A > r* and thus modifications to their growth. 

To obtain the precise value for r* in a given 
cosmology one has to solve the nonlinear equation 
Eq. (81) (which is meant to represent Eq. 59). The 
answer for the Vainshtein scale r, can only be given 
in a statistical sense, that is, the solution of that 
equation will depend on the spectrum of perturba- 
tions and cosmological parameters. In the context 
of the Schwarzschild solution, one can think of it as 
having to estimate r, from a collection of objects 
(how many and of what masses depending on the 
spectrum of perturbations) in a situation where the 
superposition principle does not apply (hence it is 
not possible to use the known answer for a single 
object). How different objects interact is controlled 
by the nonlinearities. We will tackle the calcula- 
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tion of the Vainshtcin scale using perturbation the- 
ory below, and N-body simulations in paper II [48]. 
See [90] for a recent discussion of the impact of Vain- 
shtein's mechanism for macroscopic violations of the 
equivalence principle. 



the nonlinear dynamics of the brane-bcnding mode 
into the Poisson equation and thus find an effective 
Poisson equation that relates the Newtonian poten- 
tial to the density perturbation. Equation (59) at 
the brane becomes (generalizing Eq. 70), 



C. The Effective Poisson Equation 

Let us now discuss in detail the behavior of pertur- 
bations including nonlinearities. We can incorporate 



(87) 



which we can write in Fourier space as, 



- '^^'1 + '^^ + ^ k^Co + Hqak $0 + fc' ' 



1 

2^ 



(fcl • - HH Co(fel)Co(fe2) [6d] d''kid'k2, 



where [fc] = (5D(fc — fci — ^2)- This, together with the Poisson equation, Eq. (71), which in view of Eq. (41) 
we can rewrite in terms of ^0 as 



(89) 



determines how gravity responds to density perturbations. After expressing all metric variables interms of 
(^0 and Co, using that $0 = 2^o + [v{q ~ 1) + l]C'o/'"c) we can use Eq. (89) to eliminate the brane-bending 
mode Co in terms of <i)^) and 5 and thus Eq. (88) results in an effective nonlinear Poisson equation that relates 
the gravitational potential (po to the density contrast 5, 



-(k5 + k^o) = (^) I A5 - 2(^) Y [1 - (^1 ■ ^2)'] (a^ + kf<l^o) (a5 + kl^o) [Su] d^kA 



(90) 



where 



6G = G^s-G, A = AnGa^p, 



(91) 



and we have used the fact that when the nonlinear 
term becomes important the normal derivative terms 
can be neglected, i.e. a/{krc) <C 1 at small scales. 

It's easy to see that this equation simplifies con- 
siderably in the spherical approximation. Indeed, in 
that case the kernel [1 — (fci • ^2)^] =2/3 and Fourier 



transforming back to real space, 
local quadratic equation for (V^i 
S, whose solution is 



3A/5G\ 
4g\G J 



Eq. (90) becomes a 
io — A(5) with source 



1+3.^-1 



(92) 
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where g is a dimensionless nonlinearity parameter 
given by 



~ \G J \ aJ 3(2?72 - 277 + 1)^ ' 



9 = 2 



(93) 



and decreases as time goes on (due to the dilution 
in matter density) from g = 1/3 at early times to 
g « 0.1 at z = 0. As a result, the transition to 
standard gravity happens at relatively large den- 
sity contrasts S > 10. One can also interpret the 
spherical modified Poisson equation, by defining a 
^-dependent Newton's constant by 



G 



A 



dS 



jSG/G) 



(94) 



which also emphasizes that for voids, ^ < 0, is 
further suppressed (recall SG/G < in our case) 
compared to (5 > 0, that's another way of saying 
that new non-Gaussianities will be induced. 

Equation (92) is the same as the modified Pois- 
son equation given by Eq. (3.1) in [47], after using 
that their variables e and can be written in our 
notation as e = {8/3)gS and = 3 (SG/G). Note, 
however, that our treatment includes the scale de- 
pendence of SG/G induced by the normal derivatives 
to the branc at the linear level, neglected in [47]. 

In paper II we also perform numerical simulations 
under the spherical approximation (used in [91] to 
simulate DGP and degravitation theories) and com- 
pare them to the full solution to check their validity. 
We find that the spherical approximation is accurate 
to about 10% in the density power spectrum. 

Equation (90) gives, beyond the spherical approx- 
imation, a nonlocal and nonlinear relationship be- 
tween V^^o and the density perturbations 6. We will 
explore two consequences of this equation. First, we 
discuss how to calculate the relationship between c/jq 
and S, i.e. the propagator for Eq. (90), using an ap- 
proximate resummation technique that captures the 
non-local corrections beyond the spherical approxi- 
mation. We will use this to calculate the nonlinear 
power spectrum and evaluate the transition scale to 
GR. Finally, we use perturbation theory in Eq. (90) 
to calculate the newly induced non-Gaussianities in 
the bispectrum. 



D. The Modified Poisson Propagator 

The right hand side of Eq. (90) describes the cor- 
rections to the standard Poisson equation. At large 
scales, density perturbations are small and the first, 
linear, term dominates over the second, nonlinear, 
contribution. This is what leads to the modified 
growth of perturbations described by an effective 
scale and time dependent gravitational constant. As 
small scales are approached, however, the nonlinear 
term becomes important, and as we discussed above 
it makes the right hand side of Eq. (90) subdomi- 
nant compared to the left hand side, recovering GR 
at small scales. 

To follow the nonlinear growth of structure we are 
interested in how V^^o responds to density fluctu- 
ations 5. Since we are interested in calculating sta- 
tistical averages (e.g. the power spectrum) we can 
describe the response in Fourier space as 

^5,5,^3 _ l^s,S2 ) ^3 - eye.) 



3\\VSiVS2V63 



+ 



(95) 



where V denotes a functional derivative, integration 
over repeated Fourier arguments is understood, and 
Si = 5{ki). In the language of renormalized pertur- 
bation theory [92, 93] the coefficient of the linear 
term can be regarded as the two-point propagator 
for the modified Poisson equation, 

r,(fc)fc(fe-fe')-(^^>, (96) 

the coefficient of the quadratic term as its three- 
point propagator 



rf (fci,fe2)fc(fe-fei2) 



1 / v^Mk) 



2\\VS{ki)VS{k2)/' 
(97) 

and so on. The two-point propagator describes 
the renormalized linear response, with asymptotics 



r^{kr» < 1) 



A Geff 

A;2 G ' 

r^(fcn>i) ~ 



(98) 
(99) 
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FIG. 4: Top panel: The two-point propagator for the 
modified Poisson equation in the naive spherical approx- 
imation, Eq. (100), for difi'erent redshifts corresponding 
to z = (sohd), z = 0.75 (long dashed), and z = 1.5 
(short dashed). Bottom panel: Same for the three-point 
propagator, Eq. (101). Note that the time dependence 
is not monotonic in this case. 



where is the Vainshtein scale. Recall that at large 
scales the effective Newton's constant GcS depends 
on fc, with nonlocal terms giving T(f,(k ^ 0) oc rcfc~^, 
in precise correspondence with the graviton propa- 
gator. On the other hand, at small scales where ex- 
trinsic curvature fluctuations are subdominant com- 
pared to density (intrinsic curvature) fluctuations 
returns to its GR (unmodifled Poisson) value. The 
transition happens through nonlinear effects in the 
dynamics of the brane bending mode, which is the 
nontrivial part of evaluating the expectation value 
in Eq. (96). 

(2) 

The three-point propagator is important at 
large scales where it gives rise to additional non- 
Gaussianity. However, as we shall see below it is a 
rather small correction (few percent) at large scales, 
whereas at small scales we know from the discus- 
sion on the Vainshtein mechanism above that it 



must be highly suppressed since corrections to lin- 
ear response at small scales go roughly as extrin- 
sic curvature fluctuations which scale as ^ 5^^^, 
which from Eq. (95) suggest F^^^ ~ A~'^/^, where 

A^(fc) = Ank^Pik) is the amplitude of density per- 
turbations. Indeed, it is easy to see this in the spher- 
ically symmetric case. From Eq. (92) we can roughly 
estimate 



r0(fc) 



A 



{5G/G) 



1 + |5A(fc) 



(100) 



where we naively replaced 5 — > A(fc) when taking 
expectation values in Eq. (94). This is of course not 
a valid calculation of the spherical propagator (the 
dependence on the power spectrum is simply a guess, 
and higher-order statistics than the power spectrum 
should enter into the correct spherical answer), but 
it is enough for estimation purposes. To make this 
distinction clear we refer to this calculation as the 
naive spherical approximation. Similarly, we can es- 
timate the three-point propagator as 



(fc) 



A 2 f5G\ 



l+3ffA(fc) 



-3/2 



(101) 



Thus we see that when A 1, F 



(2) 





A-3/2 and 



that F^ ^ F^ at all scales. The overall scale of 

(2) 

F^ ' compared to F^ is set by the nonlinear cou- 
pling 5, which is always much smaller than unity 
(see Eq. 93). Higher-order propagators will be sup- 
pressed by correspondingly more powers of mak- 
ing the description of the modifled Poisson equation 
by the two-point propagator alone a very good ap- 
proximation. 

Figure 4 shows these results as a function of k for 
redshifts z — 0,0.75,1, where for the power spec- 
trum we have used the HaloFit fltting formula [94] 
for standard gravity corresponding to a linear nor- 
malization given by the DGP expansion history. We 
will compare these simple results with the calcula- 
tion including effects beyond the spherical approxi- 
mation that follows below. 

All these arguments suggest that to obtain an ef- 
fective resummation of the two-point propagator F^ 
is very reasonable to use linear response in Eq. (90). 
For convenience in writing the equations that follow 



20 



below we write the propagator in terms of a new 
function M{k), 

M{k) + k'^Mk) = -A (^) (fc) 6{k) M{k), (102) 
thus 

A r 1 

(103) 



r0(fc) = -^ + 



where M describes on average the Vainshtein non- 
linear dynamics, with asymptotics 



M(fc^O)^l, M(fc^oo) 



0. 



(104) 



We must note that the resummation is essential to 
recover the large-fc limit, and that M{k) in Eq. (102) 
denotes effective (in mean field sense) corrections to 
the Poisson equation, which is what enters when cal- 
culating correlation functions; i.e. it is not meant to 



describe a particular realization in the ensemble (for 
that it would have to be a random field itself). We 
may regard the function M{k) as a form-factor that 
describes the suppression of the brane-bending mode 
at small scales; or, in terms of the gravitoii propa- 
gator, it corresponds to including such momentum 
dependent form-factor in the second term on the sec- 
ond line of Eq. (22), corresponding to the hclicity- 
zero mode, thus leading to the disappearance of the 
vDVZ discontinuity at high momentum. 

The nonlinearity of the effective Poisson equation 
is encoded in the fact that M depends on the spec- 
trum and bispcctrum of density perturbations, while 
the non-locality essentially comes from the fact that 
the bispcctrum has a nontrivial dependence on tri- 
angle shape. Although we only explicitly used k as 
an argument, M also depends on time, as do all vari- 
ables in Eqs. (102-103). In order to find an equation 
for M we rewrite Eq. (90) using Eq. (102), 



(M - 1) 6{k) 



1 - (fci • SiMi 62M2 [(5d] SkiSk2, 



(105) 



and multiplying by 5{k'), taking expectation values, and Fourier transforming back to real space, we find 



iir) - x(r) = g j 
where the correlation functions are. 



1 - (fci • hf] B{ki,k2) MiM2e'^'^-^ d^kid^k2, 



^{r)= j ^k-"" P{k)d'k, 



e'^-^ P{k)M{k) d?k, 



and P and B are the power spectrum and bispcctrum of density perturbations, 

( 5{k)5{k') ) = P{k) Sn{k + k'), { d{ki)6{k2)6{k3) ) = B{kuk2) ^(fci + fes + feg). 



(106) 



(107) 



(108) 



As mentioned above, the non-Gaussian corrections induced by the nonlincaritics in the modified gravity 
sector are small, since they are suppressed by powers of g (see Section V C below for an explicit calculation of 
the bispcctrum). Therefore, for the purpose of calculating M(fc) we shall use the standard gravity bispcctrum. 
Since the Vainshtein scale is in the nonlinear regime [47] , we need a description of the bispcctrum at nonlinear 
scales. Thus we use the bispcctrum fitting formula obtained in [95] , which reads 

B{kuk2)=2F2{ki,k2)PlP2+2F2{k2,k3)P2P3 + 2F2ik3,ki)P3Pi, (109) 

where fcs = —ki2 and the kernel F2 is given by 

F2{kuk2) = ^ a{ki)a{k2) + ^(^i • k2){^ + ^) h{kx)b{k2) + ^(fci • kf c(fci)c(fc2), (HO) 
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where the functions a, b and c arc given by [95] 



_ l+a8"° V0-7Q3N (g/4)"+3-5 1 + 0.4 (n + 3) _ 1 + 4.5/[1.5 + (n + 3)^] (2g)"+3 

l + (g/4)"+3-5 ' ~ l + gn+3.5 ' l + (2g)»+3-5 

(111) 

where g = k/kni with A;„i the scale at which the dimensionless hnear power is unity, n is the effective spectral 
index at scale k, and Qz{n) = (4 — 2")/(l + 2"+^). At large scales, a,b,c ^ 1 and tree-level perturbation 
theory is recovered, at small scales b,c ^ and a hierarchical bispcctrum with saturation value Q^in) 
follows [96] . Possible deviations from the hierarchical ansatz in the highly nonlinear regime will only impact 
M{k) at very large wavenumbers, where it is nearly zero. 

In order to solve for M{k) from Eq. (106) it is convenient to define an effective amplitude Qes{r), 

B{ki,k2) M1M2 e*'^!^-^ (fkid^k2, (112) 



Qcff(r) [x(r)] = 



1 - (fci • k2f 



which converts Eq. (106) into a local equation for 
X{r), 

ar)-x{r)=gQes{r)[x{r)?, (113) 

and the solution for M is obtained from the result- 
ing quadratic equation for x{f) by inverse Fourier 
transform, 



M{k) 



(27r 



2.gOoffP(fc) 



(114) 

It's easy to see that this M{k) satisfies the asymp- 
totics in Eq. (104) and represents an effective re- 
summation of the Vainshtein mechanism. The 
Vainshtein scale for cosmological perturbation thus 
depends on the characteristic nonlinearity g (see 
Eq. 93), the bispectrum of perturbations (since the 
nonlinearities in the field equations are quadratic) 
weighted by the nonlocal kernel of second deriva- 
tives 1 — (fci • ^2)^ represented here by Qeff, and the 
power spectrum of perturbations. 

To carry out the calculation of M, we then proceed 
as follows: 

i) We start assuming the naive spherical solution 
for the modified Poisson propagator (and thus 
M), Eq. (100). 

ii) Given this M, we calculate QeS from Eq. (112) 

at the desired redshift. For the power spec- 
trum, we use the nonlinear power spectrum 
given by HaloFit [94] corresponding to a linear 



spectrum in standard gravity with the growth 
factor that follows the modified expansion his- 
tory, let us call this GRH. For example, for the 
Vtm = 0.2 model we use, we assume trg = 0.9 at 
= for GR, then cts = 0.754 for GRH, and 
<T8 = 0.696 for linearized DGP. Since the Qeff 
ratio involves equal number of power spectra 
in numerator and denominator, the assumed 
power spectrum should not affect the answer 
by much. As we shall see, the GRH power 
spectrum never differs from the fully nonlin- 
ear DGP model by more than 30%. 

iii) Evaluate the new M{k) from Eq. (114). 

iv) With this new solution for M{k) we go back of 
step i) and iterate these steps until convergence 
is achieved. 

Figure 5 shows the results of this iterating proce- 
dure for z = 0. Short dashed lines show the initial 
guess that corresponds to the naive spherical solu- 
tion for M (bottom panel) and that implied for Q^s 
(top panel) after using Eq. (112), long dashed lines 
shows the solutions after 3 iterations, and solid lines 
after 10 iterations. We see that convergence is quite 
fast, typically a few iterations (at high redshift con- 
vergence is much faster) . Some noise develops in Q^s 
at small scales (r < lh~^ Mpc), but as we checked 
by fitting a smooth function to Qeff) this does not 
affect the behavior of M{k). 

Figure 6 shows the effective gravitational constant 
as a function of scale when the resummation of the 
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FIG. 5; Convergence of the iteration scheme to solve the 
integral equation for the modified Poisson propagator at 
z = 0. Short dashed lines show the initial guess that 
corresponds to the naive spherical solution, long dashed 
lines shows the solutions after 3 iterations, and solid lines 
after 10 iterations. The top panel corresponds to Qcft, 
see Eqs. (112-113), while the bottom panel shows the 
function M that defines the modified Poisson propaga- 
tor, see Eqs. (103) and (114). 



modified Poisson propagator is included. We see 
that nonlinearities in the modified gravity sector 
drive Gcs/G to unity at small scales, as expected 
from the Vainshtein mechanism. We now discuss 
what this renormalization of the gravitational con- 
stant implies for the nonlinear power spectrum. 
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FIG. 6: Effective gravitational constant (solid line) as a 
function of scale at z = 0, 1, 3 after resummation of non- 
linearities in the brane-bending mode (Vainshtein mech- 
anism). Dashed lines show the linear theory calculation, 
as in Fig. 2. 



stant (two-point modified Poisson propagator) and 
what this implies for the nonlinear power spectrum, 
and second by including the three-point propaga- 
tor and what it implies for the bispectrum. Let us 
change to a time variable in terms of the scale factor 
a, to 

x = Ina, (115) 
then we can rewrite conservation of stress-energy as 



V. STATISTICS 

A. Nonlinear Equations of motion 

We now consider the equations of motion for den- 
sity and velocity fields at the nonlinear level, and 
include the nonlinearities in the modified gravity sec- 
tor, first in terms of the running gravitational con- 



— - V • V = V(,5v), (116) 

where we have written the velocity field as v = — TYv, 
and using = V • v, 

c)G - 

— + (l + g)0-V2,/,o = V-[(v-V)v], (117) 
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FIG. 7: Ratio of renormalized growth factor to standard 
linear growth factor as a function of scale for z = 0, 1. 
The enhancement towards small scales reflects the onset 
of the Vainshtein mechanism, by which gravity becomes 
GR. 



where V is the dimensionless gradient, 



and we can write the modified Poisson equation as. 



0o(fc) = r^(fc)^(fc)+ y v^f{k^M)hh [5Y,]d'kid'k2. 

(119) 

Equations (116), (117) and (119) form the closed 
system of equations to be solved, and constitutes the 
starting point to compute nonlinear corrections to 
the power spectrum and bispectrum. We can write 
a second-order equation for the evolution of density 
perturbations (putting linear terms in the left hand 
side and quadratic nonlinearities in the right hand 
side). 



d^5 



+ (i + g)|^ + Pr^<^= [(l + g) + ^ V • (,5v) + V • [(v • V)v] 



Vf 8162 [fe] d?kiSk2 (120) 



where we have slightly abused notation by writing 
some terms in real space (the standard gravity non- 
linearities) and others in Fourier space (the modified 
gravity nonlinearities, which are resummed into the 
Poisson propagators P^'s). This form of the equa- 
tions of motion is useful to derive linear and second- 
order perturbative solutions. 

At this point it is worth mentioning that had we 
carried out the full perturbative expansion of the 
DGP equations (see [97]), we would end up with an 
infinite series of vertices corresponding to the per- 
turbative solution of Eq. (90) for 00 as a function 
of (5. These terms will have two effects when cal- 



culating correlation functions: they will renormalize 
the growth factor (essentially because one can dress 
any linear propagator inside a diagram by includ- 
ing such vertices: this is described by P^), and they 

can change the mode-coupling properties (described 

(2) 

by P^ and higher-order ones). While the first ef- 
fect leads to the Vainshtein mechanism through a 
dressing of the gravitational constant, the second is 
subdoniinant to the mode-coupling effects in GR, 
because the vertices are suppressed by powers of g 
as opposed to those in GR which are always of order 
unity. Thus, while we will need the resummation of 
the two-point Poisson propagator to describe nonlin- 
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ear scales, the three-point (and higher) propagators 
only affect non-Gaussianity weakly on large scales, 
and provide very small corrections to the power spec- 
trum on intermediate scales. 



B. Power Spectrum 

From Eq. (120) we see that one effect of nonlin- 
earities in the modified gravity sector is to renor- 
malize the linear growth factor (which can be ob- 
tained from this equation by using a with M — 1 
in Eq. 103). As M decays to zero at large k (see 
Fig. 5) the effective gravitational constant runs to G 
(see Fig. 6) and gravity is enhanced, becoming GR 
in the large k limit. This means that the renormal- 
ized growth factor will increase compared to the lin- 
ear growth at small scales. Figure 7 shows the result 
for the renormalized growth factor £)™norm obtained 
by solving Eq. (120) with zero right hand side, in 
terms of the linear growth factor as a function 
of scale for z = 0, 1. We see the expected enhance- 
ment; the effect is larger at lower redshift where non- 
linearities are stronger, as expected. Because the 
two-point propagator obeys the right asymptotics 
at small scales, Eq. (104), the renormalized growth 
factor becomes the linear growth factor in GR with 
the same expansion history as DGP. That means, 
for example, that parametrizing the growth factor 
in terms of a growth index 7 [98], 7 will run with 
scale and time, becoming closer to the GR value at 
small scales and late times. 

The right hand side of Eq. (120) has the famil- 
iar form of standard gravity except for the last term 
containing the three-point propagator. As discussed 
above this is a small correction (suppressed by g) and 
decays very fast into the nonlinear regime. We shall 
quantify the size of this term in detail below when we 
discuss the bispectrum; those results indicate that 
the three-point propagator will induce extra mode 
coupling power at the 1-2% level. Therefore, within 
the accuracy of our calculations here it is well jus- 
tified to neglect this contribution when computing 
the nonlinear power spectrum. 

Under this approximation, we are thus left with 
a rather remarkable situation. The whole effect of 
nonlinearities in the modified gravity sector is to 
renormalize the growth factor, or rather renormal- 
ize the linear power spectrum. To calculate the non- 
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FIG. 8: The nonlinear density power spectrum at 2 = 
(top) and z = 1 (bottom) for the DGP model (solid lines) 
and GR with the same expansion history (dashed lines) . 



linear power spectrum, we just have to evolve this 
renormalized linear spectrum by the standard non- 
linearities we commonly deal with in GR. For this 
purpose, for simplicity we use the HaloFit [94] fit- 
ting formula that relates the linear spectrum to the 
nonlinear spectrum at given time. While this fit- 
ting formula was developed for CDM models with a 
linear evolution that is scale independent (and G in- 
dependent of scale and time), note that the time 
and scale dependence of G in our case is rather 
small, its change is only appreciable over many Hub- 
ble times and orders of magnitude in spatial scales. 
The renormalized growth factor changes even more 
slowly than G, and an adiabatic adjustment of the 
linear power spectrum at each time and scale (as 
given by Fig. 7) should be a reasonable approxima- 
tion. Note in addition that, clearly, our method 
correctly recovers large and small scales: at large 
scales nonlinearities are small and we recoverer full 
linear theory (including the transition to 5D behav- 
ior in the quasistatic approximation), at small scales 
where the Vainshtein mechanism makes gravity GR, 
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FIG. 9: The ratio of Pgrh, the power spectrum for GR 
with a modified expansion history, to that in the DGP 
model for two values of Sim , highlighting the dependence 
on the somewhat different Vainshtein scales. 



our renormalizcd growth will bring the linear power 
spectrum to the amplitude it would have in GR with 
a modified expansion history. We shall compare our 
predictions to numerical simulations in detail in pa- 
per II. 

Figure 8 shows the results for the nonlinear power 
spectrum at z = 0, 1. The dashed lines correspond 
to the GRH nonlinear spectrum, i.e. a model with 
the same expansion history (or Hubble constant) but 
with standard gravity [129]. The difference in large- 
scale normalization (erg = 0.696 versus erg = 0.754 
at z = 0) reflects the difference in gravitational 
force law alone, and is stronger at lower redshift 
where the modified gravity has had more time to 



act. This difference gets enhanced at intermediate 
scales, but this is an expected outcome of the dif- 
ference in normalizations, i.e. one-loop corrections 
would be stronger for Pgrh than for Pdgp! the ef- 
fect of the running G in Pdgp is subdominant. At 
small scales the running of G in Pdgp becomes im- 
portant and Pdgp approaches Pgrh- Some of these 
features are also seen in f{R) models, see [99]. 

Figure 9 shows more clearly these effects by tak- 
ing the ratio Pgrh /Pdgp at z = for two different 
values of ilm, the standard = 0.2 that we used 
throughout this paper and Vim = 0.27 (which we use 
in paper II to compare to simulations). Both models 
would have trg = 0.9 at z = in ACDM. The lower 
r^m is, the larger is, and thus the larger the Vain- 
shtein scale is, with stronger nonlinearities in the 
modified gravity sector (i.e. larger g). Therefore, 
although the difference between these two models is 
significant (up to 5% for k <1 ft,Mpc~^), at smaller 
scales the low Q.m model approaches Pgrh some- 
what faster (the dashed line in Fig. 9 approaches 
the solid line at large k). Our results for the behav- 
ior of Pgrh / Pdgp are qualitatively similar to those 
in [97], which use perturbation theory to estimate 
the quadratic correction to the Poisson equation and 
use it in a fitting formula proposed by [100] to in- 
terpolate between Pdgp for the linearized Poisson 
equation at large scales and Pgrh at small scales 
based on the amplitude of the dimensionless power 
A(/c) = 47rfc3p(fc). 



C. Bispectrum 

Let us now discuss the new non-Gaussianities 
coming from the modified gravity sector, which 
are induced by the extrinsic curvature fluctuations. 
Since we are interested in large scales (but not as 
large so nonlocal terms become important), we can 
use the large-scale behavior of the Poisson propaga- 
tors in Eq. (120), 



A Gcff 
fc2 G ' 



A 5G 



1 - (fci • k2 



(121) 



To compute the large-scale bispectrum, we need the second-order solution to Eq. (120). Since the new 
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contributions do not add new multipoles (only changing the i — 0,2 amplitudes), we follow the standard 
approach from PT (see e.g. [3]). We write the second-order solution as, 



5(2) (fe) = 



\k2 kj 



Liik ■ k) + ^2 L2{k ■ k) S1S2 [6d] d^kid^k2, 



(122) 



where the 1?^ are arbitrary functions of time, and Lf are Legendre polynomials. The terms inside the square 
brackets reflect the structure of the right-hand side of Eq. (120) and correspond to the kernel F2 in Eq. (110) 
in the large-scale limit {a = b = c = 1). From the constraint that ( S^^^ ) = it follows that Vi = Vq + 7^2, 
and the two independent multipole amplitudes obey the following equations of motion (a dot denotes d/dx), 



Vo + {l + q)T>o 



G 



-Vn 



SG 



G V 



(123) 
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SG/rj 
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(124) 



where D+ is the linear growth factor. These equa- 
tions are straightforward to solve numerically given 
initial conditions as determined by Eq. (110) in the 
large-scale limit, i.e. Vq — 17/21 and I?2 ~ 4/21 
with Vi oc (£'+)2. From Eq. (122) one can com- 
pute the bispectrum using Eq. (109), and from it 
the reduced bispectrum Q123 = Bi2-i/ {PiP2+P2P^ + 
P3P1). Figure 10 shows the result of the ratio of Q to 
that in standard gravity for fi^ = 0.2 and triangles 
with ki =0.1/1 Mpc~^ and = 2fci as a function of 
angle 9 between the wavcvcctors, showing that there 
is a few percent enhancement for isosceles configura- 
tions, and no correction at all for squeezed triangles 
{9 = 0,tt) where the kernel in Eq. (110) vanishes. 
The Q ratio is roughly independent of scale, and it 
is most sensitive to ftm, increasing with decreasing 
flm, again as a result that lower Qm leads to a larger 
fc and thus larger g in Eq. (121). We compare these 
results against measurements in N-body simulations 
in paper II. 

Although the amplitude of this correction is small 
compared to the current statistical errors in redshift 
surveys, the amplitude is model dependent and may 
be significantly (though only parametrically) differ- 
ent in other modifications of gravity. The character- 
istic shape dependence (no correction for squeezed 
triangles, and maximal difference for isosceles tri- 
angles) should be robust for the class of models of 



massive gravity where GR is restored by derivative 
interactions as in the Vainshtein mechanism. To the 
extent that such a mechanism is universal [45], we 
are after a generic feature of modified gravity models 
that can be tested through the study of cosmologi- 
cal perturbations. In fact, this signature is generic 
to "Galilcon" models that display a shift symmetry 
(perhaps remnant of a higher-dimensional Lorcntz 
symmetry) and second-order equations of motion, 
see [34] where all possible nonlinear operators in 4D 
with these properties arc derived (for recent discus- 
sion and extensions of the Galileon model see [101- 
104]). Indeed, it is easy to see that all these opera- 
tors do posses the same properties as far as higher- 
order correlators are concerned. That is, using the 
notation of [34], a term in the equations of motion 



(125) 



gives rise as we have seen at subhorizon scales to a 
quadratic kernel 



ki ki - (fei • k2r 



(126) 



for the bispectrum that vanishes for squeezed tri- 
angles (where fei is parallel or antiparallel to ^2). 
Similarly, the next invariant 



£4 = (acf - 3{DC){d^d,cf 



2Tr(9^5,C)3, 

(127) 
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FIG. 10: The ratio of reduced bispectra Q in brane- 
induced gravity to GR for triangles with fci = 
0.1/iMpc~^ and k2 — 2ki as a function of angle be- 
tween fci and k2. 



gives rise to a cubic kernel 
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for the trispectrum that vanishes for squeezed 
trispectrum configurations (where all wavectors are 
paraUel or antiparallel to each other). It's easy to 
show that the same is true for the last invariant £^ 
as far as the five-point function in Fourier space is 
concerned. 

Finally, it is worth noting that in the normal 
branch (instead of the self-accelerated branch as we 
worked so far), the signature in the bispectrum is 
of opposite sign (suppression for isosceles triangles), 
since (SG/G) changes sign (gravity is stronger than 
in GR). This should be a generic outcome in mod- 
els where the extra scalar degree of freedom is not a 
ghost. For a calculation of the bispectrum and skew- 
ness in f{R) models see [105] and [106], respectively. 



D. Implications for BAO, Weak Lensing 

We now comment briefly on implications of our 
results for the Baryon Acoustic Oscillations (BAO) 
method and weak gravitational lensing as tools to 
study cosmic acceleration. BAO constraints have 
been apphed to the DGP model [60, 61, 107] with 
the assumption that nonlinear effects coming from 
the modified gravity sector would not affect the de- 
termination of the acoustic scale done in observa- 
tions [108]. We are now in a better position to assess 
this assumption. 

In GR nonlinearities do shift the acoustic scale 
compared to linear theory but the effect is small, 
at the percent level [109-112]. The effect can be 
understood by the contribution of a 90-degree out 
of phase oscillation by mode-coupling effects [109], 
which develops because nonlinear corrections are a 
strong function of spectral index [113] that is mod- 
ulated by the linear spectrum BAO's. The bulk of 
the acoustic scale shift is due to the power spectrum 
of the second-order density field, Eq. (122), in par- 
ticular the dipole term £ = 1 which describes the 
(v • V) transport terms in the equation of motion, 
Eq. (120). The £ = 0, 2 terms contribute to the shift 
as the mode-coupling power has cross terms of ^ = 1 
with them. Except perhaps for creation of order 
one higher £ amplitudes, there are thus two possible 
changes to the shift in a modified theory of gravity: 
a change in the £ = 0,2 amplitudes, and a change 
on the £ = 1 amplitude itself. 

As we have seen, due to the particular form of the 
nonlinearities in the DGP model (related to the shift 
symmetry mentioned above), only £ = 0,2 are sig- 
nificantly changed, but from the calculation of the 
bispectrum, we know that this change is very small, 
at the few percent level. Therefore, such modifica- 
tion induces changes in the acoustic scale shift by a 
few percent of a percent, which is negligibly small. 
The other possible source of concern is a change in 
the £ = 1 amplitude. Such terms are indeed present 
in this theory, but are highly suppressed. They arise 
through terms such as WC-W^ in the Gauss-Codazzi 
equation, Eq. (76). It's easy to see that these are of 
order {aH/k)^ times smaller (< 10^"* at the acoustic 
scale) than the GR dipole, and that's why we have 
neglected them. 

To summarize, to appreciably change the acous- 
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tic scale through nonhncar effects beyond GR, one 
needs to induce order one changes in the £ = 0,1,2 
ampHtudes of the second-order density field, or 
equivalently, of the large-scale bispectrum. This is 
not what happens in this theory, essentially because 
the Vainshtein scale is so much smaller than the non- 
linear scale from standard gravity. In any case, one 
should be able to correct for nonlinear shifts in any 
theory of gravity by simultaneously determining the 
configuration dependence of the bispectrum with the 
power spectriim (which can determine the £ = 0,1,2 
amplitudes), as pointed out in [109]. 

One caveat needs to be brought up at this point, 
however. The implicit assumption on this BAO dis- 
cussion so far is that only mode-coupling power can 
cause shifts of the linear acoustic scale, but in theo- 
ries with scale-dependent growth factors (as generic 
gravity modifications are), another possible source 
is that the density field propagator (the smoothing 
kernel that damps linear oscillations in the power 
spectrum) may contain secondary oscillations, un- 
like the case in GR where it can be shown to be 
smooth [114]. The perceptive reader may have 
noticed that in Figs. 6 and 7, (Ges/G) and thus 
^£)renorm^^^^ present Small amplitude oscillations. 
These are about 0.2% in amplitude for the power 
spectrum (a factor of two smaller compared to the 
90-degree out of phase oscillations induced by mode- 
coupling in GR [109]). At this point it is difficult for 
us to assess the importance of this, since our cal- 
culations are definitely not accurate to that level. 
Among other things, we relied on the HaloFit fitting 
formula which has not been calibrated for spectra 
with BAO's in it. A definite assessment of this ef- 
fect will have to await the development of renormal- 
ized perturbation theory for these theories, which is 
beyond the scope of this paper. 

As far as weak lensing at subhorizon scales is con- 
cerned, as first pointed out in [47], the predictions 
for the relation of the lensing potential (0o + V'o, in 
our convention) to density perturbations S in this 
theory are indistinguishable from GR (see Eq. 43), 
as the brane-bending mode C does not couple to 
photons since it is sourced by the trace of the en- 
ergy momentum tensor and photons are conformal 
sources (equation of state w = 1/3) and, in addition, 
G exactly vanishes on the background solution, see 
Eq. (1). Therefore, weak lensing in this type of the- 
ories is completely unaffected by the extra dynamics 



in the modified gravity sector, and has no extra in- 
formation than the density field itself; i.e. the effec- 
tive Newton's constant driving the lensing potential, 
Eq. (43), is not modified. This is not strictly true 
at the largest scales, where other degrees of freedom 
enter (through the normal derivatives to the brane), 
gravity becomes 5D and this is of course followed by 
the lensing potential [130], see e.g. [2, 3, 91], this 
may be again described by a running of the gravita- 
tional constant appearing in Eq. (43), which ignored 
normal derivatives. 



VI. CONCLUSIONS 

In this paper, we have studied linear and non- 
linear solutions for cosmological perturbations in 
brane-induced gravity in five dimensions. Our re- 
sults provide improved understanding of the linear 
solutions and derive new predictions at the non- 
linear level that should be also relevant for other 
theories of modified gravity with similar character- 
istics. Most of the predictions developed here are 
compared to numerical simulations in paper II, find- 
ing good agreement. Given the length of this paper, 
it seems appropriate here to summarize the main 
results, discuss their validity and suggest possible 
improvements. 

We solved for the linear evolution of perturba- 
tions taking advantage of the symmetry under gauge 
transformations along the extra-dimension to decou- 
ple the bulk equations in the quasistatic approxima- 
tion. This is a different route than the popular ap- 
proach in the literature of using the master equation, 
and provides in our view a more physically transpar- 
ent procedure (to the extent that five- dimensional 
physics can be so!). 

We have cheeked the validity of the quasistatic 
approximation for all metric perturbations, and 
pointed out that the derivatives normal to the brane 
belong to the quasistatic approximation, giving rise 
to 5D behavior in the static on-brane potentials, 
in close analogy to the graviton propagator. In 
the past, these contributions have been considered 
as deviations from the quasistatic approximation, 
e.g. [64, 65]; therefore, the validity of this approx- 
imation for cosmological perturbations is likely bet- 
ter than thought before. In addition, it is worth 
emphasizing that these normal derivatives induce. 
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from the 4D point of view, unusual nonlocal opera- 
tors (see Eq. 69), this cannot be obtained by a local 
4D theory no matter how complicated, thus making 
brane-induced gravity distinct from arbitrary (but 
local) dark energy models. A similar regime is ex- 
pected in theories with more extra dimensions, as 
gravity "cascades down" from higher to lower di- 
mensionality as smaller distances are probed. 

Linear theory can only be, however, a limited de- 
scription of the physics in the modified gravity sec- 
tor, because if valid at all scales the theory would 
be ruled out by local experiments, e.g. by show- 
ing that light deflection and nonrelativistic motion 
display different gravitational constants. Nonlinear 
effects, through the Vainshtein mechanism, provide 
a natural way out of this situation. We studied the 
nonlinearities in the bulk and brane equations, and 
justified the approximations made, which form the 
basis for carrying out numerical simulations of struc- 
ture formation (see paper II). Along the way, we 
justified why the ghost-like brane-bending mode in- 
teracting with healthy modes does not cause catas- 
trophic instabilities for cosmological solutions. 

We studied in detail the Vainshtein mechanism, 
by which the theory becomes GR at small scales. 
In this regard, it is useful to think geometrically, in 
terms of the extrinsic curvature K of the brane; af- 
ter all, cosmic acceleration is driven in this theory by 
a balance between intrinsic (Ricci, R) and extrinsic 
curvature. At the level of perturbations, and within 
the quasistatic approximation, there are three dis- 
tinct regimes: 

1. 5D regime, with 5K/K ~ 5R/R and 5K dom- 
inated by normal derivatives, giving rise to 5D 
behavior and a scale-dependent renormaliza- 

tion of G, 

2. scalar-tensor regime, with 5K/K 5R/R and 
6K dominated by the brane bending mode C, 
leading to scale-independent (but time depen- 
dent) renormalization of the gravitational con- 
stant, 

3. Vainshtein regime, with SK/K <C SR/R and 
SK dominated by the brane bending mode C, 
but subdominant compared to density fluctu- 
ations {5K/K ~ S^/^), thus its effect on the 
physics at the brane is greatly reduced, Ges/G 



becomes scale-dependent through nonlinear ef- 
fects, approaching unity at small scales where 
GR is recovered. 

These three stages can be seen in Fig. 6, although 
the second stage is greatly reduced at low redshift. 

We showed that throughout evolution, the weak 
field approximation is always obeyed (this is also 
supported by numerical simulations in paper II) , the 
quantity that goes nonlinear in the Vainshtein mech- 
anism is the extrinsic curvature fluctuation SK/K. 
We also showed that the condition SK/K — 1 can 
be written in the familiar way (r^ = rgr^) for the 
Vainshtein scale of isolated perturbations, making 
the connection with the standard Schwarzschild so- 
lution. 

We calculated the impact of the Vainshtein mech- 
anism on the nonlinear density power spectrum. To 
do so, we derived an explicit modified Poisson equa- 
tion that relates the Newtonian potential (po to den- 
sity perturbations 6, Eq. (90). We showed that this 
nonlinear integral equation can be rewritten as an 
effective linear plus quadratic response for (j)o as 
a function of 6 (Eq. 119) by introducing the two- 
and three-point propagators of the Poisson equa- 
tion. The two-point propagator describes the run- 
ning of the gravitational constant due to the Vain- 
shtein mechanism, and in order to recover the right 
asymptotics at small scales (see Eq. 99) one has 
to resum it. The three-point propagator only con- 
tributes at large scales, and gives rise to new sig- 
natures of modified gravity as corrections to non- 
Gaussianity. 

To resum the two-point modified Poisson propa- 
gator one has to solve a complicated integral equa- 
tion, which involves the density correlators we are 
interested in; we implemented an iterative proce- 
dure based on knowledge of the power spectrum and 
bispectrum in GR with the same expansion history. 
Our method can be improved by recomputing the 
power spectrum and bispectrum at each step given 
the previous iteration solution of the two-point prop- 
agator, for self-consistency. Another aspect that 
may merit improvement is the development of bet- 
ter GR fitting functions for the power spectrum and 
bispectrum. To put our resummation in a more rig- 
oroTis footing, one needs to develop in full e.g. renor- 
malized perturbation theory [92] for modified grav- 
ity theories (also needed for checking on additional 
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BAO systomatics), which is left for future work. 

After resummation of the two-point propagator is 
achieved, the effective equations of motion for den- 
sity and velocity fields is remarkably similar to stan- 
dard gravity, see Eq. (120). The difference is a slowly 
varying gravitational constant, and a new vertex 
that gives rise to non-Gaussianity. Neglecting the 
latter (justified by the small corrections induced in 
the bispcctrum) we showed that to calculate the non- 
linear power spectrum, one can implement the GR 
fitting formula for the power spectrum on a renor- 
malized linear spectrum that already contains the 
Vainshtein mechanism through the renormalization 
of G. This renormalized linear spectrum can also be 
used for predicting the mass function of dark matter 
halos (sec paper II). The calculated nonlinear power 
spectrum (Fig. 8) approaches that of GR with the 
same expansion history, as expected, but the transi- 
tion is rather gentle. See paper II for comparison of 
these predictions with N-body simulations. 

It is worth noting that at the level of the power 
spectrum, and to a reasonably good approximation, 
the effect of nonlinearities in the modified gravity 
sector may be absorbed into a renormalization of the 
gravitational constant. Since the relation between 
the lensing potential and density perturbations is en- 
tirely unaffected by the extra physics in these theo- 
ries, the modified gravity can be described in this ap- 
proximation by a single function, an effective grav- 
itational constant Gcii{k,a) for non-rclativistic mo- 
tion that depends on space and time; the caveat is 
that it cannot be computed from linear theory. This 
simple result (checked against simulations in pa- 
per II) may have important practical consequences 
for parametrizing modified gravity in observations. 

The extra nonlinear vertex in Eq. (120) gives rise 
to additional non-Gaussianity, a signature of the 
Vainshtein mechanism at large-scales. We calcu- 
lated the corrections induced by it in the bispectrum. 



and while the amplitude is small and model depen- 
dent, the characteristic shape dependence (vanishing 
for squeezed triangles) is generic to modified gravity 
theories that are of "Galileon" type [34], i.e. where 
the dynamics of the additional scalar degree of free- 
dom is invariant under shifts of its gradient. We 
pointed out that similar results hold for higher-order 
correlators than the bispectrum in these type of the- 
ories. These signatures of nonlinearities in modified 
gravity theories can be used as a diagnostic in future 
observations. 
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